home *** CD-ROM | disk | FTP | other *** search
/ Power Programmierung / Power-Programmierung (Tewi)(1994).iso / magazine / drdobbs / 1991 / 10 / eigencal.asc < prev    next >
Text File  |  1991-09-10  |  9KB  |  310 lines
  1. _MIXED-LANGUAGE WINDOWS PROGRAMMING_
  2. by John Norwood
  3.  
  4.  
  5.  
  6. [LISTING ONE]
  7.  
  8.  
  9.       subroutine tred2(a,n,np,d,e)
  10. c
  11. c Householder reduction of a real symmetric, nxn matrix a, stored in an
  12. c npxnp physical array.  On output, a is replaced by the orthogonal matrix
  13. c q effecting the transformation.
  14. c d returns the diagonal elements of the tridiagonal matrix, and e the off-
  15. c diagonal elements, with e(1)=0.
  16.  
  17.       implicit none
  18.  
  19.       integer*4 i,j,k,l      ! Loop indices
  20.       integer*4 n            ! Actual array size
  21.       integer*4 np           ! Physical array size of incoming array
  22.       real*4    a(np,np)     ! Matrix, np is used so dimensions match
  23.       real*4    d(np)        ! Vector that will have diagonal elements
  24.       real*4    e(np)        ! Vector that will have off-diagonal elements
  25.       real*4    h            ! Vector norm used in forming projection
  26.       real*4    scale        ! Scale factor
  27.       real*4    f            ! Temporary variable
  28.       real*4    g            ! Temporary variable
  29.       real*4    hh           ! Another piece of the projection
  30.  
  31.       if (n.gt.1) then
  32.        do 18 i=n,2,-1
  33.          l=i-1
  34.          h=0.
  35.          scale=0.
  36.          if(l.gt.1) then
  37.            do 11, k=1,l
  38.              scale=scale+abs(a(i,k))
  39. 11         continue
  40.            if(scale.eq.0.) then       ! Skip transformation
  41.              e(i)=a(i,l)
  42.            else
  43.              do 12, k=1,l
  44.                a(i,k)=a(i,k)/scale    ! Use scaled a's in transformation
  45.                h=h+a(i,k)**2          ! for eigenvectors
  46. 12           continue
  47.              f=a(i,l)
  48.              g=-sign(sqrt(h),f)
  49.              e(i)=scale*g
  50.              h=h-(f*g)
  51.              a(i,l)=f-g
  52.              f=0.
  53.              do 15,j=1,l
  54.                a(j,i)=a(i,j)/h
  55.                g=0.
  56.                do 13, k=1,j
  57.                  g=g+a(j,k)*a(i,k)
  58. 13             continue
  59.                if(l.gt.j) then
  60.                do 14,k=j+1,l
  61.                  g=g+a(k,j)*a(i,k)
  62. 14             continue
  63.                endif
  64.                e(j)=g/h               ! Form element of projection in
  65.                f=f+e(j)*a(i,j)        ! temporarily unused element of e
  66. 15           continue
  67.              hh=f/(h+h)
  68.              do 17, j=1,l
  69.                f=a(i,j)
  70.                g=e(j)-hh*f
  71.                e(j)=g
  72.                do 16, k=1,j                    ! Loop to reduce matrix a
  73.                  a(j,k)=a(j,k)-f*e(k)-g*a(i,k)
  74. 16             continue
  75. 17           continue
  76.              endif
  77.            else
  78.              e(i)=a(i,l)
  79.            endif
  80.            d(i)=h
  81. 18       continue
  82.  
  83. c This starts the eigenvector specific part of the code.
  84.  
  85.        endif
  86.        d(1)=0.
  87.        e(1)=0.
  88.        do 23, i=1,n          ! Begin accumulation of transformation matrices
  89.          l=i-1
  90.          if(d(i).ne.0.) then
  91.            do 21, j=1,l
  92.              g=0.
  93.              do 19, k=1,l          ! Use information stored in a to form
  94.                g=g+a(i,k)*a(k,j)   ! projection times orthogonal matrix
  95. 19           continue
  96.              do 20, k=1,l
  97.                a(k,j)=a(k,j)-g*a(k,i)
  98. 20           continue
  99. 21         continue
  100.          endif
  101.  
  102. c This ends the eigenvector specific part of the code.
  103.  
  104.          d(i)=a(i,i)
  105.  
  106. c This starts the eigenvector specific part of the code.
  107.  
  108.          a(i,i)=1.            ! Reset row and column of a to identity matrix
  109.          if(l.ge.1)then       ! for the next iteration of transformation loop
  110.            do 22,j=1,l
  111.              a(i,j)=0.
  112.              a(j,i)=0.
  113. 22         continue
  114.          endif
  115.  
  116. c This ends the eigenvector specific part of the code.
  117.  
  118. 23     continue
  119.        return
  120.        end
  121.  
  122.       subroutine tqli(d,e,n,np,z,iterflag)
  123.  
  124. c This subroutine performs a QL algorithm with implicit shifts, to determine
  125. c the eigenvalues and eigenvectors of a real, symmetric, tridiagonal matrix,
  126. c or of a real, symmetric matrix previously reduced by tred2 above.
  127. c d is a vector of length np.  On input, its first n elements are the
  128. c diagonal elements of the tridiagonal matrix.  On output, it returns the
  129. c eigenvalues.  The vector e inputs the subdiagonal elements of the
  130. c tridiagonal matrix, with e(1) arbitrary. On output, e is destroyed.
  131. c If eigenvectors are desired, the matrix z (nxn stored in an npxnp array)
  132. c is input as the identity matrix or the matrix that is returned from tred2.
  133. c On output, the kth column of z contains the normalized eigenvector
  134. c corresponding to the eigenvalue in d(k).
  135.  
  136.       implicit none
  137.  
  138.       integer*4 i,k,l       ! Loop indices
  139.       integer*4 iter        ! Iteration counter
  140.       integer*4 n           ! Logical size of array z
  141.       integer*4 m           ! Submatrix size
  142.       integer*4 np          ! Physical size of array z
  143.       integer*4 iterflag    ! Flag to return error code to main routine
  144.       real*4    d(np)       ! Diagonal elements is, eigenvalues out
  145.       real*4    e(np)       ! Off diagonal elements in, nothing out
  146.       real*4    z(np,np)    ! Matrix from tred2 in, eigenvectors out
  147.       real*4    dd          ! Holds small subdiagonal element
  148.       real*4    g           ! Holds Givens rotation
  149.       real*4    s           ! "Sin" component of Givens rotation
  150.       real*4    c           ! "Cos" component of Givens rotation
  151.       real*4    p           ! Element of projection matrix
  152.       real*4    f           ! Holds result of "sin" applied to element of e
  153.       real*4    b           ! Holds result of "cos" applied to element of e
  154.       real*4    r           ! Temporary piece of c or s
  155.  
  156.  
  157.       if(n.gt.1) then
  158.         do 11, i=2,n                 ! Renumber elements of e
  159.           e(i-1)=e(i)
  160. 11      continue
  161.         e(n)=0.
  162.         do 15,l=1,n
  163.           iter=0
  164. 1         do 12,m=l,n-1              ! Search for small subdiagonal element
  165.             dd=abs(d(m))+abs(d(m+1))
  166.             if (abs(e(m))+dd.eq.dd) go to 2
  167. 12        continue
  168.           m=n
  169. 2         if(m.ne.l) then
  170.  
  171.             if(iter.eq.30) then
  172.              iterflag = -1
  173.              return
  174.             endif
  175.  
  176.             iter=iter+1
  177.             g=(d(l+1)-d(l))/(2.*e(l))   ! Calculate shift
  178.             r=sqrt(g**2+1.)
  179.             g=d(m)-d(l)+e(l)/(g+sign(r,g))
  180.             s=1.
  181.             c=1.
  182.             p=0.
  183.             do 14, i=m-1,l,-1           ! Plane rotation followed by a Givens
  184.               f=s*e(i)                  ! rotations to maintain tridiagonal
  185.               b=c*e(i)                  ! form
  186.               if(abs(f).ge.abs(g))then
  187.                 c=g/f
  188.                 r=sqrt(c**2+1.)
  189.                 e(i+1)=f*r
  190.                 s=1./r
  191.                 c=c*s
  192.               else
  193.                 s=f/g
  194.                 r=sqrt(s**2+1.)
  195.                 e(i+1)=g*r
  196.                 c=1./r
  197.                 s=s*c
  198.               endif
  199.               g=d(i+1)-p
  200.               r=(d(i)-g)*s+2.*c*b
  201.               p=s*r
  202.               d(i+1)=g+p
  203.               g=c*r-b
  204.  
  205. c Start of code specific to forming eigenvectors.
  206.  
  207.               do 13, k=1,n
  208.                 f=z(k,i+1)
  209.                 z(k,i+1)=s*z(k,i)+c*f
  210.                 z(k,i)=c*z(k,i)-s*f
  211. 13            continue
  212.  
  213. c End of code specific to forming eigenvectors.
  214.  
  215. 14          continue
  216.             d(l)=d(l)-p
  217.             e(l)=g
  218.             e(m)=0.
  219.             go to 1
  220.           endif
  221. 15        continue
  222.         endif
  223.  
  224.  
  225.         return
  226.         end
  227.  
  228.  
  229.  
  230.  
  231. Examplσ 1:
  232.  
  233. (a)
  234.  
  235.        SUBROUTINE STRINGER(INSTRING)
  236.        CHARACTER*40 INSTRING
  237.        INSTRING = 'This is from Fortran'
  238.        END
  239.  
  240. (b)
  241. è
  242.  
  243. Sub Form_Click ()
  244.   Dim temp As String * 40
  245.   Call STRINGER(temp)
  246.   Debug.Print temp
  247. End Sub
  248.  
  249. (c)
  250.  
  251. Declarσ SuΓ STRINGE╥ LiΓ "d:\vb\test\string.dlló (ByVa∞ Mystrinτ A≤ String⌐ 
  252.  
  253.  
  254.  
  255.  
  256. Examplσ 2:
  257.  
  258. (a)
  259.        SUBROUTINE ARRAYTEST(ARR)
  260.        INTEGER*4 ARR(20)
  261.        ARR = 5
  262.        END
  263.  
  264. (b)
  265.  
  266. Sub Form_Click ()
  267.   Static testarray(1 To 20) As Long
  268.   Call ARRAYTEST(testarray(1))
  269.   For i% = 1 To 20
  270.     Debug.Print testarray(i%)
  271.   Next i%
  272. End Sub
  273.  
  274.  
  275.  
  276. (c)
  277. Declarσ SuΓ ARRAYTES╘ LiΓ "d:\vb\test\array.dlló (Myarra∙ A≤ Long)
  278.  
  279.  
  280.  
  281.  
  282. Examplσ 3
  283.  
  284. (a)
  285.  
  286.       SUBROUTINE ARRAYSTRINGS(ARR)
  287.       CHARACTER*24 ARR(5)
  288.       ARR = 'This is a string also'
  289.       END
  290.  
  291. (b)
  292.  
  293. Sub Form_Click ()
  294.   Static testarray(1 To 5) As StringArray
  295.   Call ARRAYTEST(testarray(1))
  296.   For i% = 1 To 5
  297.     Debug.Print testarray(i%).strings
  298.   Next i%
  299. End Sub
  300.  
  301.  
  302.  
  303. (c)
  304.  
  305. Type StringArray
  306.    strings As String * 24
  307. End Type
  308.  
  309. Declare Sub ARRAYSTRINGS Lib "d:\vb\test\array.dll" (Myarray As StringArray)
  310.